###Justin H Kirkland      ####
###Jonathan B Slapin      ####
###Strategic Party Disloyalty#
##In the US House         ####
###11-6-13                ####
##############################

##Clear memory
rm(list=ls())

#Load library; set seed
library(RCurl)
library(foreign)
library(arm)

#Load data
party.unity<-read.csv("House_Party_Unity_35-112.csv")
party.unity<-subset(party.unity, party.unity$Congress>=95)
votetype<-read.dta("House for Crespin Duke.dta")
vote95<-read.csv("House Roll Call 95.csv")
for(i in 10:length(names(vote95))){
  names(vote95)[i]<-paste("V", as.numeric(strsplit(names(vote95)[i], split="var")[[1]][2])-9, sep="")
}
vote96<-read.csv("House Roll Call 96.csv")
vote97<-read.csv("House Roll Call 97.csv")
vote98<-read.csv("House Roll Call 98.csv")
vote99<-read.csv("House Roll Call 99.csv")
vote100<-read.csv("House Roll Call 100.csv")
vote101<-read.csv("House Roll Call 101.csv")
vote102<-read.csv("House Roll Call 102.csv")
vote103<-read.csv("House Roll Call 103.csv")
vote104<-read.csv("House Roll Call 104.csv")
vote105<-read.csv("House Roll Call 105.csv")
vote106<-read.csv("House Roll Call 106.csv")
for(i in 8:length(names(vote106))){
  names(vote106)[i]<-paste("V", as.numeric(strsplit(names(vote106)[i], split="var")[[1]][2])-9, sep="")
}
vote107<-read.csv("House Roll Call 107.csv")
for(i in 10:length(names(vote107))){
  names(vote107)[i]<-paste("V", as.numeric(strsplit(names(vote107)[i], split="var")[[1]][2])-9, sep="")
}
vote108<-read.csv("House Roll Call 108.csv")
vote109<-read.csv("House Roll Call 109.csv")
vote110<-read.csv("House Roll Call 110.csv")
#Correct Names in 110
for(i in 10:length(names(vote110))){
  names(vote110)[i]<-paste("V", as.numeric(strsplit(names(vote110)[i], split="var")[[1]][2])-9, sep="")
}
vote111<-read.csv("House Roll Call 111.csv")
vote112<-read.csv("House Roll Call 112.csv")


#subset vote type data to 109 and party unity votes
vote.type<-subset(votetype, votetype$cong>=95 & votetype$v16==1)

#subset further to procedural and final passage votes
procedural.votes<-subset(vote.type, vote.type$vote==61 | vote.type$vote==76 |vote.type$vote==77 | vote.type$vote==79 | vote.type$vote==81 | vote.type$vote==82 | vote.type$vote==88 | vote.type$vote==89 | vote.type$vote==93 | vote.type$vote==95 | vote.type$vote==97 | vote.type$vote==99)
mispro.votes<-subset(vote.type, vote.type$vote==56 | vote.type$vote==58 | vote.type$vote==59 | vote.type$vote==66 | vote.type$vote==67 | vote.type$vote==69 | vote.type$vote==73 | vote.type$vote==74 | vote.type$vote==75 | vote.type$vote==83 | vote.type$vote==84 | vote.type$vote==86 | vote.type$vote==87 | vote.type$vote==90 | vote.type$vote==91 | vote.type$vote==92 | vote.type$vote==96 | vote.type$vote==98)
final.passage<-subset(vote.type, vote.type$vote==11 | vote.type$vote==12 | vote.type$vote==14 | vote.type$vote==30)

#places to store data
party.unity$PartyUnityPro<-numeric(nrow(party.unity))
party.unity$PartyVotesPro<-numeric(nrow(party.unity))
party.unity$PartySupportPro<-numeric(nrow(party.unity))

party.unity$PartyUnityMisPro<-numeric(nrow(party.unity))
party.unity$PartyVotesMisPro<-numeric(nrow(party.unity))
party.unity$PartySupportMisPro<-numeric(nrow(party.unity))

party.unity$PartyUnityPass<-numeric(nrow(party.unity))
party.unity$PartyVotesPass<-numeric(nrow(party.unity))
party.unity$PartySupportPass<-numeric(nrow(party.unity))


#loop to fill in data
for(i in 1:nrow((party.unity))){
  #for(i in 1:15){   #for testing
  year<-party.unity$Congress[i]   #ID congress of observation
  
  #Relevent roll call votes
  if(year==95){votes<-vote95}
  if(year==96){votes<-vote96}
  if(year==97){votes<-vote97}
  if(year==98){votes<-vote98}
  if(year==99){votes<-vote99}
  if(year==100){votes<-vote100}
  if(year==101){votes<-vote101}
  if(year==102){votes<-vote102}
  if(year==103){votes<-vote103}
  if(year==104){votes<-vote104}
  if(year==105){votes<-vote105}
  if(year==106){votes<-vote106}
  if(year==107){votes<-vote107}
  if(year==108){votes<-vote108}
  if(year==109){votes<-vote109}
  if(year==110){votes<-vote110}
  if(year==111){votes<-vote111}
  if(year==112){votes<-vote112}
  
  #New matrix with just voting
  vote.mat<-votes[,10:ncol(votes)]
  
  #Match observation to roll call row
  w<-which(votes$id==party.unity$ICPSRid[i])
  
  #Locate pro and pass votes in roll call matrix
  provotes<-procedural.votes$voteview[which(procedural.votes$cong==year)]
  passvotes<-final.passage$voteview[which(final.passage$cong==year)]
  mprovotes<-mispro.votes$voteview[which(mispro.votes$cong==year)]
  
  #Correct names
  names(vote.mat)<-gsub("V", "", names(vote.mat))
  
  ###################
  #Conditional on Party Affiliation
  if(party.unity$party[i]==100 | party.unity$party[i]==200){
    
    
    ###############################################
    #Calculate Party Support on Final Passage######
    party.unity$PartyVotesPass[i]<-sum(vote.mat[w,match(passvotes, names(vote.mat))]!=9) #Party Line Votes
    
    #Find Party choices on party line votes
    pv<-apply(vote.mat[which(votes$party==party.unity$party[i]), passvotes], 2, table)
    pv2<-numeric(length(pv))
    for(j in 1:length(pv)){
      pv2[j]<-as.numeric(names(which.max(pv[[j]])))  
    }
    
    #How many times did observation support party
    party.unity$PartySupportPass[i]<-sum(vote.mat[w,match(passvotes, names(vote.mat))]==pv2 & pv2!=9)
    
    #Calculate percentage
    party.unity$PartyUnityPass[i]<-(party.unity$PartySupportPass[i]/party.unity$PartyVotesPass[i])*100
    
    #####################################################
    #Calculate Party Support on Partisan Procedure#######
    party.unity$PartyVotesPro[i]<-sum(vote.mat[w,match(provotes, names(vote.mat))]!=9) #Party Line Votes
    
    #Find Party choices on party line votes
    pv<-apply(vote.mat[which(votes$party==party.unity$party[i]), provotes], 2, table)
    pv2<-numeric(length(pv))
    for(j in 1:length(pv)){
      pv2[j]<-as.numeric(names(which.max(pv[[j]])))  
    }
    
    #How many times did observation support party
    party.unity$PartySupportPro[i]<-sum(vote.mat[w,match(provotes, names(vote.mat))]==pv2 & pv2!=9)
    
    #Calculate percentage
    party.unity$PartyUnityPro[i]<-(party.unity$PartySupportPro[i]/party.unity$PartyVotesPro[i])*100
    
    #####################################################
    #Calculate Party Support on Misc Procedure#######
    party.unity$PartyVotesMisPro[i]<-sum(vote.mat[w,match(mprovotes, names(vote.mat))]!=9) #Party Line Votes
    
    #Find Party choices on party line votes
    pv<-apply(vote.mat[which(votes$party==party.unity$party[i]), mprovotes], 2, table)
    pv2<-numeric(length(pv))
    for(j in 1:length(pv)){
      pv2[j]<-as.numeric(names(which.max(pv[[j]])))  
    }
    
    #How many times did observation support party
    party.unity$PartySupportMisPro[i]<-sum(vote.mat[w,match(mprovotes, names(vote.mat))]==pv2 & pv2!=9)
    
    #Calculate percentage
    party.unity$PartyUnityMisPro[i]<-(party.unity$PartySupportMisPro[i]/party.unity$PartyVotesMisPro[i])*100
  }
  
  
  ###################
  #Conditional on Party Affiliation
  if(party.unity$party[i]!=100 & party.unity$party[i]!=200){
    party.unity$PartyVotesPass[i]<-NA
    party.unity$PartySupportPass[i]<-NA
    party.unity$PartyUnityPass[i]<-NA
    party.unity$PartyVotesPro[i]<-NA
    party.unity$PartySupportPro[i]<-NA
    party.unity$PartyUnityPro[i]<-NA
    party.unity$PartyVotesMisPro[i]<-NA
    party.unity$PartySupportMisPro[i]<-NA
    party.unity$PartyUnityMisPro[i]<-NA
  }
  print(i)  #for tracking
}


######################################
#Read in Nominate Data################
nom.data<-read.dta("NomScores.DTA")

##############################################
####Merge Full Data###########################
nom.data$Congress<-nom.data$cong
nom.data$ICPSRid<-nom.data$idno

master2<-merge(party.unity, nom.data, by=c("ICPSRid", "Congress"))

########Lag Nominate Scores
names(master2)
master2$lag.nom<-rep(NA, nrow(master2))
for(i in 1:nrow(master2)){
  z<-which(master2$ICPSRid==master2$ICPSRid[i] & master2$Congress==(master2$Congress[i]-1))
  master2$lag.nom[i]<-ifelse(length(z)>0, master2$dwnom1[z], NA)
  print(i)
}


####################################################
###Interactions with Majority Party Status##########
master.d<-subset(master2, master2$party.x==100)
master.r<-subset(master2, master2$party.x==200)

#ID majority party years
master.d$majority<-ifelse(master.d$Congress==95|master.d$Congress==96|master.d$Congress==97|master.d$Congress==98|master.d$Congress==99|master.d$Congress==100|master.d$Congress==101|master.d$Congress==102|master.d$Congress==103|master.d$Congress==110|master.d$Congress==111, 1, 0)
master.r$majority<-ifelse(master.r$Congress==104|master.r$Congress==105|master.r$Congress==106|master.r$Congress==107|master.r$Congress==108|master.r$Congress==109|master.r$Congress==112, 1, 0)

#Identify Congresses
Congress<-sort(unique(master.d$Congress))

################################################
###Difference of means test#####################
master.d$Ex<-numeric(nrow(master.d))
master.d$Mod<-numeric(nrow(master.d))
for(i in 1:length(Congress)){
  z<-quantile(abs(master.d$dwnom1[which(master.d$Congress==Congress[i])]), 0.9) 
  z2<-quantile(abs(master.d$dwnom1[which(master.d$Congress==Congress[i])]), 0.1) 
  q<-which(master.d$Congress==Congress[i])  
  master.d$Ex[q]<-z
  master.d$Mod[q]<-z2
}

master.d$Extremist<-ifelse(abs(master.d$dwnom1)>=master.d$Ex,1, NA)
master.d$Extremist<-ifelse(abs(master.d$dwnom1)<=master.d$Mod,0, master.d$Extremist)

#####TABLE 1: Dems##########
###Diff for Extremist vs Moderates
t.test(master.d$PartyUnityPass[which(master.d$Extremist==1)]~master.d$majority[which(master.d$Extremist==1)])
t.test(master.d$PartyUnityPass[which(master.d$Extremist==0)]~master.d$majority[which(master.d$Extremist==0)])

master.r$Ex<-numeric(nrow(master.r))
master.r$Mod<-numeric(nrow(master.r))
for(i in 1:length(Congress)){
  z<-quantile(abs(master.r$dwnom1[which(master.r$Congress==Congress[i])]), 0.9) 
  z2<-quantile(abs(master.r$dwnom1[which(master.r$Congress==Congress[i])]), 0.1) 
  q<-which(master.r$Congress==Congress[i])  
  master.r$Ex[q]<-z
  master.r$Mod[q]<-z2
}

master.r$Extremist<-ifelse(abs(master.r$dwnom1)>=master.r$Ex,1, NA)
master.r$Extremist<-ifelse(abs(master.r$dwnom1)<=master.r$Mod,0, master.r$Extremist)

#####TABLE 1: Reps##########
###Diff for Extremist vs Moderates
t.test(master.r$PartyUnityPass[which(master.r$Extremist==1)]~master.r$majority[which(master.r$Extremist==1)])
t.test(master.r$PartyUnityPass[which(master.r$Extremist==0)]~master.r$majority[which(master.r$Extremist==0)])

########################################
################Table 2#################
#Fractional Logit models
summary(mod<-glmer(matrix(c(PartySupport, PartyVotes-PartySupport), ncol=2)~abs(lag.nom)+majority+majority:abs(lag.nom)+(1|Congress)+(1|ICPSRid), family=binomial(link=logit), data=master.d, control=glmerControl(optimizer="bobyqa")))
summary(mod2<-glmer(matrix(c(PartySupport, PartyVotes-PartySupport), ncol=2)~abs(lag.nom)+majority+majority:abs(lag.nom)+(1|Congress)+(1|ICPSRid), family=binomial(link=logit), data=master.r, control=glmerControl(optimizer="bobyqa")))

#Marginal effects
MFX.d1<-fixef(mod)[2]+fixef(mod)[4]*c(0,1)
se.MFX.d1<-1.96*sqrt(vcov(mod)[2,2]+c(0,1)^2*vcov(mod)[4,4]+2*c(0,1)*vcov(mod)[2,4])
MFX.r1<-fixef(mod2)[2]+fixef(mod2)[4]*c(0,1)
se.MFX.r1<-1.96*sqrt(vcov(mod2)[2,2]+c(0,1)^2*vcov(mod2)[4,4]+2*c(0,1)*vcov(mod2)[2,4])

#####################################
#Sims for Figures####################

library(mvtnorm)
sims<-1000
dist<-rmvnorm(sims, fixef(mod), as.matrix(vcov(mod)))
dist2<-rmvnorm(sims, fixef(mod2), as.matrix(vcov(mod2)))
calc<-300
z<-seq(min(abs(master.d$lag.nom), na.rm=T), max(abs(master.d$lag.nom), na.rm=T), length=calc)
z2<-seq(min(abs(master.r$lag.nom), na.rm=T), max(abs(master.r$lag.nom), na.rm=T), length=calc)
lo1<-numeric(calc)
hi1<-numeric(calc)
pe1<-numeric(calc)
lo2<-numeric(calc)
hi2<-numeric(calc)
pe2<-numeric(calc)
lo3<-numeric(calc)
hi3<-numeric(calc)
pe3<-numeric(calc)
lo4<-numeric(calc)
hi4<-numeric(calc)
pe4<-numeric(calc)
for(i in 1:calc){
  pp<-numeric(sims)
  pp2<-numeric(sims)
  pp3<-numeric(sims)
  pp4<-numeric(sims)
  for(j in 1:sims){
    #Dems; Min
    pp[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*0+dist[j,4]*0*z[i]
    #Dems; Maj
    pp2[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*1+dist[j,4]*1*z[i]
    #Reps; Min
    pp3[j]<-dist2[j,1]+dist2[j,2]*z2[i]+dist2[j,3]*0+dist2[j,4]*0*z2[i]
    #Reps; Maj
    pp4[j]<-dist2[j,1]+dist2[j,2]*z2[i]+dist2[j,3]*1+dist2[j,4]*1*z2[i]
    
  }
  
  pe1[i]<-invlogit(mean(pp))
  lo1[i]<-invlogit(quantile(pp, 0.025))
  hi1[i]<-invlogit(quantile(pp, 0.975))
  
  pe2[i]<-invlogit(mean(pp2))
  lo2[i]<-invlogit(quantile(pp2, 0.025))
  hi2[i]<-invlogit(quantile(pp2, 0.975))
  
  pe3[i]<-invlogit(mean(pp3))
  lo3[i]<-invlogit(quantile(pp3, 0.025))
  hi3[i]<-invlogit(quantile(pp3, 0.975))
  
  pe4[i]<-invlogit(mean(pp4))
  lo4[i]<-invlogit(quantile(pp4, 0.025))
  hi4[i]<-invlogit(quantile(pp4, 0.975))
  
  print(i)
}



#########################
#Figure 2

#pdf("PredProbAllVotes.pdf", width=10)
par(mfrow=c(1,2))
plot(c(0.2, 0.7), c(pe1[1], pe3[1]), xlim=c(0, 1), ylim=c(0.4, 1), pch=16, col=c("blue", "red"), ylab="Predicted Probability of Party Loyalty", xlab="Party Affiliation", xaxt="n", bty="n", main="Predicted Probability of Party \nLoyalty Vote for Moderate Partisans on \nAll Roll Call Votes", cex=1.5)
grid()
points(c(0.3, 0.8), c(pe2[1], pe4[1]), pch=15, col=c("blue", "red"), cex=1.5)
segments(0.2, lo1[1], 0.2, hi1[1], lty=1, lwd=3, col="blue")
segments(0.3, lo2[1], 0.3, hi2[1], lty=1, lwd=3, col="blue")
segments(0.7, lo3[1], 0.7, hi3[1], lty=1, lwd=3, col="red")
segments(0.8, lo4[1], 0.8, hi4[1], lty=1, lwd=3, col="red")
legend("bottomleft", inset=0.05, c("Democrats", "Republicans", "Minority Party", "Majority Party"), col=c("blue", "red", "black", "black"), lty=c(1,1,NA, NA), pch=c(NA, NA, 16, 15), cex=0.8, lwd=c(3,3,NA,NA))

plot(c(0.2, 0.7), c(pe1[300], pe3[300]), xlim=c(0, 1), ylim=c(0.97, 1), pch=16, col=c("blue", "red"), ylab="Predicted Probability of Party Loyalty", xlab="Party Affiliation", xaxt="n", bty="n", main="Predicted Probability of Party \nLoyalty Vote for Extreme Partisans on \nAll Roll Call Votes", cex=1.5)
points(c(0.3, 0.8), c(pe2[300], pe4[300]), pch=15, col=c("blue", "red"), cex=1.5)
grid()
segments(0.2, lo1[300], 0.2, hi1[300], lty=1, lwd=3, col="blue")
segments(0.3, lo2[300], 0.3, hi2[300], lty=1, lwd=3, col="blue")
segments(0.7, lo3[300], 0.7, hi3[300], lty=1, lwd=3, col="red")
segments(0.8, lo4[300], 0.8, hi4[300], lty=1, lwd=3, col="red")
#dev.off()

##########################
########add in pro's######
#Table 3 Part 1###########
summary(mod<-glmer(matrix(c(PartySupportPro, PartyVotesPro-PartySupportPro), ncol=2)~abs(lag.nom)+majority+majority:abs(lag.nom)+(1|Congress)+(1|ICPSRid), family=binomial(link=logit), data=master.d, control=glmerControl(optimizer="bobyqa")))
summary(mod2<-glmer(matrix(c(PartySupportPro, PartyVotesPro-PartySupportPro), ncol=2)~abs(lag.nom)+majority+majority:abs(lag.nom)+(1|Congress)+(1|ICPSRid), family=binomial(link=logit), data=master.r, control=glmerControl(optimizer="bobyqa")))

#Marginal Effects
MFX.d3<-fixef(mod)[2]+fixef(mod)[4]*c(0,1)
se.MFX.d3<-1.96*sqrt(vcov(mod)[2,2]+c(0,1)^2*vcov(mod)[4,4]+2*c(0,1)*vcov(mod)[2,4])
MFX.r3<-fixef(mod2)[2]+fixef(mod2)[4]*c(0,1)
se.MFX.r3<-1.96*sqrt(vcov(mod2)[2,2]+c(0,1)^2*vcov(mod2)[4,4]+2*c(0,1)*vcov(mod2)[2,4])


library(mvtnorm)
sims<-1000
dist<-rmvnorm(sims, fixef(mod), as.matrix(vcov(mod)))
dist2<-rmvnorm(sims, fixef(mod2), as.matrix(vcov(mod2)))
calc<-300
z<-seq(min(abs(master.d$lag.nom), na.rm=T), max(abs(master.d$lag.nom), na.rm=T), length=calc)
z2<-seq(min(abs(master.r$lag.nom), na.rm=T), max(abs(master.r$lag.nom), na.rm=T), length=calc)
lo1<-numeric(calc)
hi1<-numeric(calc)
pe1<-numeric(calc)
lo2<-numeric(calc)
hi2<-numeric(calc)
pe2<-numeric(calc)
lo3<-numeric(calc)
hi3<-numeric(calc)
pe3<-numeric(calc)
lo4<-numeric(calc)
hi4<-numeric(calc)
pe4<-numeric(calc)
for(i in 1:calc){
  pp<-numeric(sims)
  pp2<-numeric(sims)
  pp3<-numeric(sims)
  pp4<-numeric(sims)
  for(j in 1:sims){
    #Dems; Min
    pp[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*0+dist[j,4]*0*z[i]
    #Dems; Maj
    pp2[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*1+dist[j,4]*1*z[i]
    #Reps; Min
    pp3[j]<-dist2[j,1]+dist2[j,2]*z2[i]+dist2[j,3]*0+dist2[j,4]*0*z2[i]
    #Reps; Maj
    pp4[j]<-dist2[j,1]+dist2[j,2]*z2[i]+dist2[j,3]*1+dist2[j,4]*1*z2[i]
    
  }
  
  pe1[i]<-invlogit(mean(pp))
  lo1[i]<-invlogit(quantile(pp, 0.025))
  hi1[i]<-invlogit(quantile(pp, 0.975))
  
  pe2[i]<-invlogit(mean(pp2))
  lo2[i]<-invlogit(quantile(pp2, 0.025))
  hi2[i]<-invlogit(quantile(pp2, 0.975))
  
  pe3[i]<-invlogit(mean(pp3))
  lo3[i]<-invlogit(quantile(pp3, 0.025))
  hi3[i]<-invlogit(quantile(pp3, 0.975))
  
  pe4[i]<-invlogit(mean(pp4))
  lo4[i]<-invlogit(quantile(pp4, 0.025))
  hi4[i]<-invlogit(quantile(pp4, 0.975))
  
  print(i)
}


#############################
#Figure 3 Part 1

#pdf("PredProbProVotes.pdf", width=10)
par(mfrow=c(1,2))
plot(c(0.2, 0.7), c(pe1[1], pe3[1]), xlim=c(0, 1), ylim=c(0.4, 1), pch=16, col=c("blue", "red"), ylab="Predicted Probability of Party Loyalty", xlab="Party Affiliation", xaxt="n", bty="n", main="Predicted Probability of Party \nLoyalty for Moderate Partisans on \nProcedural Roll Call Votes", cex=1.5)
grid()
points(c(0.3, 0.8), c(pe2[1], pe4[1]), pch=15, col=c("blue", "red"), cex=1.5)
segments(0.2, lo1[1], 0.2, hi1[1], lty=1, lwd=3, col="blue")
segments(0.3, lo2[1], 0.3, hi2[1], lty=1, lwd=3, col="blue")
segments(0.7, lo3[1], 0.7, hi3[1], lty=1, lwd=3, col="red")
segments(0.8, lo4[1], 0.8, hi4[1], lty=1, lwd=3, col="red")
legend("bottomleft", inset=0.05, c("Democrats", "Republicans", "Minority Party", "Majority Party"), col=c("blue", "red", "black", "black"), lty=c(1,1,NA, NA), pch=c(NA, NA, 16, 15), cex=0.8, lwd=c(3,3,NA,NA))


plot(c(0.2, 0.7), c(pe1[300], pe3[300]), xlim=c(0, 1), ylim=c(0.97, 1), pch=16, col=c("blue", "red"), ylab="Predicted Probability of Party Loyalty", xlab="Party Affiliation", xaxt="n", bty="n", main="Predicted Probability of Party \nLoyalty for Extreme Partisans on \nProcedural Roll Call Votes", cex=1.5)
points(c(0.3, 0.8), c(pe2[300], pe4[300]), pch=15, col=c("blue", "red"), cex=1.5)
grid()
segments(0.2, lo1[300], 0.2, hi1[300], lty=1, lwd=3, col="blue")
segments(0.3, lo2[300], 0.3, hi2[300], lty=1, lwd=3, col="blue")
segments(0.7, lo3[300], 0.7, hi3[300], lty=1, lwd=3, col="red")
segments(0.8, lo4[300], 0.8, hi4[300], lty=1, lwd=3, col="red")
#dev.off()

############################
######Final Passage#########

#Table 3 Part 2#############
summary(mod<-glmer(matrix(c(PartySupportPass, PartyVotesPass-PartySupportPass), ncol=2)~abs(lag.nom)+majority+majority:abs(lag.nom)+(1|Congress)+(1|ICPSRid), family=binomial(link=logit), data=master.d, control=glmerControl(optimizer="bobyqa")))
summary(mod2<-glmer(matrix(c(PartySupportPass, PartyVotesPass-PartySupportPass), ncol=2)~abs(lag.nom)+majority+majority:abs(lag.nom)+(1|Congress)+(1|ICPSRid), family=binomial(link=logit), data=master.r, control=glmerControl(optimizer="bobyqa")))

#Marginal Effects
MFX.d2<-fixef(mod)[2]+fixef(mod)[4]*c(0,1)
se.MFX.d2<-1.96*sqrt(vcov(mod)[2,2]+c(0,1)^2*vcov(mod)[4,4]+2*c(0,1)*vcov(mod)[2,4])
MFX.r2<-fixef(mod2)[2]+fixef(mod2)[4]*c(0,1)
se.MFX.r2<-1.96*sqrt(vcov(mod2)[2,2]+c(0,1)^2*vcov(mod2)[4,4]+2*c(0,1)*vcov(mod2)[2,4])


library(mvtnorm)
sims<-1000
dist<-rmvnorm(sims, fixef(mod), as.matrix(vcov(mod)))
dist2<-rmvnorm(sims, fixef(mod2), as.matrix(vcov(mod2)))
calc<-300
z<-seq(min(abs(master.d$lag.nom), na.rm=T), max(abs(master.d$lag.nom), na.rm=T), length=calc)
z2<-seq(min(abs(master.r$lag.nom), na.rm=T), max(abs(master.r$lag.nom), na.rm=T), length=calc)
lo1<-numeric(calc)
hi1<-numeric(calc)
pe1<-numeric(calc)
lo2<-numeric(calc)
hi2<-numeric(calc)
pe2<-numeric(calc)
lo3<-numeric(calc)
hi3<-numeric(calc)
pe3<-numeric(calc)
lo4<-numeric(calc)
hi4<-numeric(calc)
pe4<-numeric(calc)
for(i in 1:calc){
  pp<-numeric(sims)
  pp2<-numeric(sims)
  pp3<-numeric(sims)
  pp4<-numeric(sims)
  for(j in 1:sims){
    #Dems; Min
    pp[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*0+dist[j,4]*0*z[i]
    #Dems; Maj
    pp2[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*1+dist[j,4]*1*z[i]
    #Reps; Min
    pp3[j]<-dist2[j,1]+dist2[j,2]*z2[i]+dist2[j,3]*0+dist2[j,4]*0*z2[i]
    #Reps; Maj
    pp4[j]<-dist2[j,1]+dist2[j,2]*z2[i]+dist2[j,3]*1+dist2[j,4]*1*z2[i]
    
  }
  
  pe1[i]<-invlogit(mean(pp))
  lo1[i]<-invlogit(quantile(pp, 0.025))
  hi1[i]<-invlogit(quantile(pp, 0.975))
  
  pe2[i]<-invlogit(mean(pp2))
  lo2[i]<-invlogit(quantile(pp2, 0.025))
  hi2[i]<-invlogit(quantile(pp2, 0.975))
  
  pe3[i]<-invlogit(mean(pp3))
  lo3[i]<-invlogit(quantile(pp3, 0.025))
  hi3[i]<-invlogit(quantile(pp3, 0.975))
  
  pe4[i]<-invlogit(mean(pp4))
  lo4[i]<-invlogit(quantile(pp4, 0.025))
  hi4[i]<-invlogit(quantile(pp4, 0.975))
  
  print(i)
}


#############################
#Figure 3 Part 2#############

#pdf("PredProbFinalVotes.pdf", width=10)
par(mfrow=c(1,2))
plot(c(0.2, 0.7), c(pe1[1], pe3[1]), xlim=c(0, 1), ylim=c(0, 1), pch=16, col=c("blue", "red"), ylab="Predicted Probability of Party Loyalty", xlab="Party Affiliation", xaxt="n", bty="n", main="Predicted Probability of Party \nLoyalty for Moderate Partisans on \nFinal Passage Roll Call Votes", cex=1.5)
grid()
points(c(0.3, 0.8), c(pe2[1], pe4[1]), pch=15, col=c("blue", "red"), cex=1.5)
segments(0.2, lo1[1], 0.2, hi1[1], lty=1, lwd=3, col="blue")
segments(0.3, lo2[1], 0.3, hi2[1], lty=1, lwd=3, col="blue")
segments(0.7, lo3[1], 0.7, hi3[1], lty=1, lwd=3, col="red")
segments(0.8, lo4[1], 0.8, hi4[1], lty=1, lwd=3, col="red")
legend("bottomleft", inset=0.05, c("Democrats", "Republicans", "Minority Party", "Majority Party"), col=c("blue", "red", "black", "black"), lty=c(1,1,NA, NA), pch=c(NA, NA, 16, 15), cex=0.8, lwd=c(3,3,NA,NA))


plot(c(0.2, 0.7), c(pe1[300], pe3[300]), xlim=c(0, 1), ylim=c(0.94, 1), pch=16, col=c("blue", "red"), ylab="Predicted Probability of Party Loyalty", xlab="Party Affiliation", xaxt="n", bty="n", main="Predicted Probability of Party \nLoyalty for Extreme Partisans on \nFinal Passage Roll Call Votes", cex=1.5)
points(c(0.3, 0.8), c(pe2[300], pe4[300]), pch=15, col=c("blue", "red"), cex=1.5)
grid()
segments(0.2, lo1[300], 0.2, hi1[300], lty=1, lwd=3, col="blue")
segments(0.3, lo2[300], 0.3, hi2[300], lty=1, lwd=3, col="blue")
segments(0.7, lo3[300], 0.7, hi3[300], lty=1, lwd=3, col="red")
segments(0.8, lo4[300], 0.8, hi4[300], lty=1, lwd=3, col="red")
#dev.off()

#############################################
###Electoral Costs of Party Loyalty##########
###Replications in R#########################
#############################################

#Clear memory
rm(list=ls())

#libraries
library(arm)
library(foreign)
library(plm)

#Load data
data<-read.dta("CKL data 1956-2004 7-23-07.dta")

#Generate some interactions
data$dwunity<-data$dwnom1recode*data$partyunity_1
data$majunity<-data$partyunity_1*data$majparty
data$dwmaj<-data$dwnom1recode*data$majparty
data$dwmajunity<-data$dwnom1recode*data$partyunity_1*data$majparty

#Gen state -dist id
data$statedist<-paste(data$icpsrst, data$cd)


#set data as plm
data.t<-pdata.frame(data, index=c("statedist", "congress"))

##############################
#Table 4#####################
#run analysis###############
summary(mod<-plm(inc2sh3_100~lag(inc2sh3_100)+partyunity_1+majparty+dwnom1recode+dwunity+majunity+dwmaj+dwmajunity+presvote+chql+spendgap+frosh+presapp_rescale_in+midterm_in+rdi_q3_in+inparty, data=data.t, effect="individual"))

#MFX
print(min.mod.mfx<-coef(mod)[2])
print(min.ex.mfx<-coef(mod)[2]+coef(mod)[5])
print(maj.mod.mfx<-coef(mod)[2]+coef(mod)[6])
print(maj.ex.mfx<-coef(mod)[2]+coef(mod)[5]+coef(mod)[6]+coef(mod)[8])

#############################
##Predicted Vote Share
library(mvtnorm)
sims<-100
calc<-150
pu<-seq(min(data$partyunity_1, na.rm=T), max(data$partyunity_1, na.rm=T), length=calc)
#pu<-seq(0.75, 1.0, length=calc)
sim.coef<-rmvnorm(sims, as.numeric(coef(mod)), as.matrix(vcov(mod)))
lo<-numeric(calc)
hi<-numeric(calc)
pe<-numeric(calc)
lo2<-numeric(calc)
hi2<-numeric(calc)
pe2<-numeric(calc)
lo3<-numeric(calc)
hi3<-numeric(calc)
pe3<-numeric(calc)
lo4<-numeric(calc)
hi4<-numeric(calc)
pe4<-numeric(calc)
for(i in 1:calc){
  pp<-numeric(sims)
  pp2<-numeric(sims)
  pp3<-numeric(sims)
  pp4<-numeric(sims)
  for(j in 1:sims){
    #Min; Mod
    pp[j]<-sim.coef[j,1]*mean(lag(data$inc2sh3_100), na.rm=T)+sim.coef[j,2]*pu[i]+sim.coef[j,3]*0+sim.coef[j,4]*0+sim.coef[j,5]*0*pu[i]+sim.coef[j,6]*0*pu[i]+sim.coef[j,7]*0*0+sim.coef[j,8]*0*0*pu[i]+sim.coef[j,9]*mean(data$presvote, na.rm=T)+sim.coef[j,10]*mean(data$chql, na.rm=T)+sim.coef[j,11]*mean(data$spendgap, na.rm=T)+sim.coef[j,12]*mean(data$frosh, na.rm=T)+sim.coef[j,13]*mean(data$presapp_rescale_in, na.rm=T)+sim.coef[j,14]*mean(data$midterm_in, na.rm=T)+sim.coef[j,15]*mean(data$rdi_q3_in, na.rm=T)+sim.coef[j,16]*mean(data$inparty, na.rm=T)
    #Min; Ex
    pp2[j]<-sim.coef[j,1]*mean(lag(data$inc2sh3_100), na.rm=T)+sim.coef[j,2]*pu[i]+sim.coef[j,3]*0+sim.coef[j,4]*1+sim.coef[j,5]*1*pu[i]+sim.coef[j,6]*0*pu[i]+sim.coef[j,7]*0*1+sim.coef[j,8]*0*1*pu[i]+sim.coef[j,9]*mean(data$presvote, na.rm=T)+sim.coef[j,10]*mean(data$chql, na.rm=T)+sim.coef[j,11]*mean(data$spendgap, na.rm=T)+sim.coef[j,12]*mean(data$frosh, na.rm=T)+sim.coef[j,13]*mean(data$presapp_rescale_in, na.rm=T)+sim.coef[j,14]*mean(data$midterm_in, na.rm=T)+sim.coef[j,15]*mean(data$rdi_q3_in, na.rm=T)+sim.coef[j,16]*mean(data$inparty, na.rm=T)
    #Maj; Mod
    pp3[j]<-sim.coef[j,1]*mean(lag(data$inc2sh3_100), na.rm=T)+sim.coef[j,2]*pu[i]+sim.coef[j,3]*1+sim.coef[j,4]*0+sim.coef[j,5]*0*pu[i]+sim.coef[j,6]*1*pu[i]+sim.coef[j,7]*0*1+sim.coef[j,8]*0*1*pu[i]+sim.coef[j,9]*mean(data$presvote, na.rm=T)+sim.coef[j,10]*mean(data$chql, na.rm=T)+sim.coef[j,11]*mean(data$spendgap, na.rm=T)+sim.coef[j,12]*mean(data$frosh, na.rm=T)+sim.coef[j,13]*mean(data$presapp_rescale_in, na.rm=T)+sim.coef[j,14]*mean(data$midterm_in, na.rm=T)+sim.coef[j,15]*mean(data$rdi_q3_in, na.rm=T)+sim.coef[j,16]*mean(data$inparty, na.rm=T)
    #Maj; Ex
    pp4[j]<-sim.coef[j,1]*mean(lag(data$inc2sh3_100), na.rm=T)+sim.coef[j,2]*pu[i]+sim.coef[j,3]*1+sim.coef[j,4]*1+sim.coef[j,5]*1*pu[i]+sim.coef[j,6]*1*pu[i]+sim.coef[j,7]*1*1+sim.coef[j,8]*1*1*pu[i]+sim.coef[j,9]*mean(data$presvote, na.rm=T)+sim.coef[j,10]*mean(data$chql, na.rm=T)+sim.coef[j,11]*mean(data$spendgap, na.rm=T)+sim.coef[j,12]*mean(data$frosh, na.rm=T)+sim.coef[j,13]*mean(data$presapp_rescale_in, na.rm=T)+sim.coef[j,14]*mean(data$midterm_in, na.rm=T)+sim.coef[j,15]*mean(data$rdi_q3_in, na.rm=T)+sim.coef[j,16]*mean(data$inparty, na.rm=T)
    #print(j)
  }
  pe[i]<-mean(fixef(mod))+mean(pp)
  lo[i]<-mean(fixef(mod))+quantile(pp, 0.025)
  hi[i]<-mean(fixef(mod))+quantile(pp, 0.975)
  pe2[i]<-mean(fixef(mod))+mean(pp2)
  lo2[i]<-mean(fixef(mod))+quantile(pp2, 0.025)
  hi2[i]<-mean(fixef(mod))+quantile(pp2, 0.975)
  pe3[i]<-mean(fixef(mod))+mean(pp3)
  lo3[i]<-mean(fixef(mod))+quantile(pp3, 0.025)
  hi3[i]<-mean(fixef(mod))+quantile(pp3, 0.975)
  pe4[i]<-mean(fixef(mod))+mean(pp4)
  lo4[i]<-mean(fixef(mod))+quantile(pp4, 0.025)
  hi4[i]<-mean(fixef(mod))+quantile(pp4, 0.975)
  print(i)
}


#####################
#Figure 4############

#pdf("VoteShareExtremist2.pdf", width=12)
par(mfrow=c(1,2))
plot(pu, pe2, ylim=c(25,75), type="l", lwd=4, ylab="Predicted Share of the Two Party Vote", xlab="Party Loyalty", bty="n")
lines(pu, lo2, lty=2, lwd=3)
lines(pu, hi2, lty=2, lwd=3)
grid()
abline(h=50)
legend("bottomright", inset=0.05, c("Point Estimate", "95% Confidence Interval"), col=c("black", "black"), lwd=c(4,3), lty=c(1,2), cex=0.8)
title(main="Ideologically Extreme Incumbents' Share \nof Two-Party Vote as Party Loyalty \nChanges (Minority Party Members)")
text(0.25, 30, "*45% Share of the Two-Party Vote", cex=0.75)
text(0.8, 70, "*64% Share of the Two-Party Vote", cex=0.75)
rug(data$partyunity_1)
plot(pu, pe4, lty=1, col="black", ylim=c(40, 80), type="l", lwd=4, ylab="Predicted Share of the Two Party Vote", xlab="Party Loyalty", bty="n")
lines(pu, lo4, lty=2, col="black", lwd=3)
lines(pu, hi4, lty=2, col="black", lwd=3)
grid()
abline(h=50)
title(main="Ideologically Extreme Incumbents' Share \nof Two-Party Vote as Party Loyalty \nChanges (Majority Party Members)")
text(0.25, 60, "*70% Share of the Two-Party Vote", cex=0.75)
text(0.8, 75, "*65% Share of the Two-Party Vote", cex=0.75)
rug(data$partyunity_1)
#dev.off()
